{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Physical modeling\n",
    "\n",
    "Code examples from [Think Complexity, 2nd edition](http://greenteapress.com/wp/complexity2), Chapter 7\n",
    "\n",
    "Copyright 2016 Allen Downey, [MIT License](http://opensource.org/licenses/MIT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function, division\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import thinkplot\n",
    "\n",
    "from matplotlib import rc\n",
    "rc('animation', html='html5')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from thinkstats2 import RandomSeed\n",
    "RandomSeed(17)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Diffusion\n",
    "\n",
    "Before we get to a Reaction-Diffusion model, we'll start with simple diffusion.\n",
    "\n",
    "The kernel computes the difference between each cell and the sum of its neighbors.\n",
    "\n",
    "At each time step, we compute this difference, multiply by a constant, and add it back in to the array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from Cell2D import Cell2D, Cell2DViewer\n",
    "from scipy.signal import correlate2d\n",
    "\n",
    "class Diffusion(Cell2D):\n",
    "    \"\"\"Diffusion Cellular Automaton.\"\"\"\n",
    "\n",
    "    kernel = np.array([[0, 1, 0],\n",
    "                       [1,-4, 1],\n",
    "                       [0, 1, 0]])\n",
    "\n",
    "    def __init__(self, n, m=None, r=0.1):\n",
    "        \"\"\"Initializes the attributes.\n",
    "\n",
    "        n: number of rows\n",
    "        m: number of columns\n",
    "        r: diffusion rate constant\n",
    "        \"\"\"\n",
    "        self.r = r\n",
    "        m = n if m is None else m\n",
    "        self.array = np.zeros((n, m), np.float)\n",
    "\n",
    "    def step(self):\n",
    "        \"\"\"Executes one time step.\"\"\"\n",
    "        c = correlate2d(self.array, self.kernel, mode='same')\n",
    "        self.array += self.r * c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We don't really have to customize the viewer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class DiffusionViewer(Cell2DViewer):\n",
    "    cmap = plt.get_cmap('Reds')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a simple example starting with an \"island\" of material in the middle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diff = Diffusion(n=9, r=0.1)\n",
    "diff.add_cells(3, 3, '111', '111', '111')\n",
    "viewer = DiffusionViewer(diff)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the initial condition:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "anim = viewer.animate()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's how it behaves over time: the \"material\" spreads out until the level is equal on the whole array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "anim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diff = Diffusion(n=9, r=0.1)\n",
    "diff.add_cells(3, 3, '111', '111', '111')\n",
    "viewer = DiffusionViewer(diff)\n",
    "\n",
    "thinkplot.preplot(cols=3)\n",
    "viewer.draw()\n",
    "\n",
    "thinkplot.subplot(2)\n",
    "viewer.step(5)\n",
    "viewer.draw()\n",
    "\n",
    "thinkplot.subplot(3)\n",
    "viewer.step(10)\n",
    "viewer.draw()\n",
    "\n",
    "plt.savefig('chap07-1.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reaction-Diffusion\n",
    "\n",
    "Now we'll add a second material and let them interact.\n",
    "\n",
    "The following function helps with setting up the initial conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def island(a, val, noise=None):\n",
    "    \"\"\"Adds an island in the middle of the array.\n",
    "            \n",
    "    val: height of the island\n",
    "    noise: magnitude of random noise\n",
    "    \"\"\"\n",
    "    noise = val if noise is None else noise\n",
    "    n, m = a.shape\n",
    "    r = min(n, m) // 20\n",
    "    a[n//2-r:n//2+r, m//2-r:m//2+r] = val\n",
    "    a += noise * np.random.random((n, m))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the RD model, we have two arrays, one for each chemical.\n",
    "\n",
    "Following [Sims](http://www.karlsims.com/rd.html), I'm using a kernel that includes the diagonal elements.  They have lower weights because they are farther from the center cell.\n",
    "\n",
    "The `step` function computes these functions:\n",
    "\n",
    "$\\Delta A = r_a \\nabla^2 A - AB^2 + f (1-A) $\n",
    "\n",
    "$\\Delta B = r_b \\nabla^2 B + AB^2 - (k+f) B $\n",
    "\n",
    "where $\\nabla^2$ is the Laplace operator the kernel is intended to approximate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class ReactionDiffusion(Cell2D):\n",
    "    \"\"\"Reaction-Diffusion Cellular Automaton.\"\"\"\n",
    "\n",
    "    kernel = np.array([[.05, .2, .05],\n",
    "                       [ .2, -1, .2],\n",
    "                       [.05, .2, .05]])\n",
    "\n",
    "    options = dict(mode='same', boundary='wrap')\n",
    "\n",
    "    def __init__(self, n, m=None, params=(0.5, 0.25, 0.035, 0.057)):\n",
    "        \"\"\"Initializes the attributes.\n",
    "\n",
    "        n: number of rows\n",
    "        m: number of columns\n",
    "        params: tuple of (Da, Db, f, k)\n",
    "        \"\"\"\n",
    "        \n",
    "        self.params = params\n",
    "        m = n if m is None else m\n",
    "        self.array = np.ones((n, m), dtype=float)\n",
    "\n",
    "        self.array2 = np.zeros((n, m), dtype=float)\n",
    "        island(self.array2, val=0.1, noise=0.1)\n",
    "        \n",
    "    def step(self):\n",
    "        \"\"\"Executes one time step.\"\"\"\n",
    "        A = self.array\n",
    "        B = self.array2\n",
    "        ra, rb, f, k = self.params\n",
    "        \n",
    "        cA = correlate2d(A, self.kernel, **self.options)\n",
    "        cB = correlate2d(B, self.kernel, **self.options)\n",
    "        reaction = A * B**2\n",
    "        self.array += ra * cA - reaction + f * (1-A) \n",
    "        self.array2 += rb * cB + reaction - (f+k) * B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The viewer for the CA shows both arrays with some transparency, so we can see where one, the other, or both, levels are high.\n",
    "\n",
    "Unlike previous CAs, the state of each cell is meant to represent a continuous quantity, so it is appropriate to interpolate.\n",
    "\n",
    "Note that `draw` has to make copies of the arrays because `step` updates the arrays in place."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class RDViewer(Cell2DViewer):\n",
    "    \"\"\"Generates images and animations.\"\"\"\n",
    "    \n",
    "    cmapu = plt.get_cmap('Reds')\n",
    "    cmapv = plt.get_cmap('Blues')\n",
    "\n",
    "    options = dict(alpha=0.7,\n",
    "                  interpolation='bicubic')\n",
    "        \n",
    "    def __init__(self, viewee):\n",
    "        \"\"\"Initializes the attributes.\n",
    "        \n",
    "        viewee: the object to be represented\n",
    "        \"\"\"\n",
    "        self.viewee = viewee\n",
    "        self.imu = None\n",
    "        self.imv = None\n",
    "        self.hlines = None\n",
    "        self.vlines = None\n",
    "\n",
    "    def draw(self, grid=False):\n",
    "        \"\"\"Draws the cells.\"\"\"\n",
    "        au = self.viewee.array.copy()\n",
    "        av = self.viewee.array2.copy()\n",
    "        \n",
    "        n, m = av.shape\n",
    "        plt.axis([0, m, 0, n])\n",
    "        plt.xticks([])\n",
    "        plt.yticks([])\n",
    "\n",
    "        self.options['extent'] = [0, m, 0, n]\n",
    "        self.imu = plt.imshow(au, cmap=self.cmapu, **self.options)\n",
    "        self.imv = plt.imshow(av, cmap=self.cmapv, **self.options)\n",
    "\n",
    "    def animate_func(self, i):\n",
    "        \"\"\"Draws one frame of the animation.\"\"\"\n",
    "        if i > 0:\n",
    "            self.step(iters=100)\n",
    "\n",
    "        self.imu.set_array(self.viewee.array)\n",
    "        self.imv.set_array(self.viewee.array2)\n",
    "        return (self.imu, self.imv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's an example using `params3`, which yields blue dots that seem to undergo mitosis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "params1 = 0.5, 0.25, 0.035, 0.057   # white spots\n",
    "params2 = 0.5, 0.25, 0.055, 0.062   # coral\n",
    "params3 = 0.5, 0.25, 0.039, 0.065   # blue spots\n",
    "\n",
    "rd = ReactionDiffusion(100, params=params3)\n",
    "viewer = RDViewer(rd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a random starting condition with lots of A, a sprinkling of B everywhere, and an island of B in the middle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "anim = viewer.animate(frames=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's how it behaves over time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "anim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I'll use the following function to generate figures using different parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_viewer(f, k, n=100):\n",
    "    \"\"\"Makes a ReactionDiffusion viewer with given parameters.\n",
    "    \"\"\"\n",
    "    params = 0.5, 0.25, f, k\n",
    "    rd = ReactionDiffusion(n, params=params)\n",
    "    viewer = RDViewer(rd)\n",
    "    return viewer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following parameters yield pink stripes and spots on a blue background:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "viewer = make_viewer(0.035, 0.057)\n",
    "\n",
    "thinkplot.preplot(cols=3)\n",
    "viewer.step(1000)\n",
    "viewer.draw()\n",
    "\n",
    "thinkplot.subplot(2)\n",
    "viewer.step(2000)\n",
    "viewer.draw()\n",
    "\n",
    "thinkplot.subplot(3)\n",
    "viewer.step(4000)\n",
    "viewer.draw()\n",
    "\n",
    "plt.savefig('chap07-2.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following parameters yield blue stripes on a pink background."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "viewer = make_viewer(0.055, 0.062)\n",
    "\n",
    "thinkplot.preplot(cols=3)\n",
    "viewer.step(1000)\n",
    "viewer.draw()\n",
    "\n",
    "thinkplot.subplot(2)\n",
    "viewer.step(2000)\n",
    "viewer.draw()\n",
    "\n",
    "thinkplot.subplot(3)\n",
    "viewer.step(4000)\n",
    "viewer.draw()\n",
    "\n",
    "plt.savefig('chap07-3.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following parameters yield blue dots on a pink background"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "viewer = make_viewer(0.039, 0.065)\n",
    "\n",
    "thinkplot.preplot(cols=3)\n",
    "viewer.step(1000)\n",
    "viewer.draw()\n",
    "\n",
    "thinkplot.subplot(2)\n",
    "viewer.step(2000)\n",
    "viewer.draw()\n",
    "\n",
    "thinkplot.subplot(3)\n",
    "viewer.step(4000)\n",
    "viewer.draw()\n",
    "\n",
    "plt.savefig('chap07-4.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Percolation\n",
    "\n",
    "In the percolation model, each cell is porous with probability `p`.  We start with a row of wet cells at the time.  During each time step, a cell becomes wet if it is porous and at least one neighbor is wet (using a 4-cell neighborhood).  For each value of `p` we compute the probability that water reaches the bottom row.\n",
    "\n",
    "Porous cells have state `1` and wet cells have state `5`, so if a cell has a wet neighbor, the sum of the neighbors will by `5` or more.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.signal import correlate2d\n",
    "\n",
    "class Percolation(Cell2D):\n",
    "    \"\"\"Percolation Cellular Automaton.\"\"\"\n",
    "\n",
    "    kernel = np.array([[0, 1, 0],\n",
    "                       [1, 0, 1],\n",
    "                       [0, 1, 0]])\n",
    "\n",
    "    def __init__(self, n, m=None, p=0.5, seed=None):\n",
    "        \"\"\"Initializes the attributes.\n",
    "\n",
    "        n: number of rows\n",
    "        m: number of columns\n",
    "        p: probability of porousness\n",
    "        \"\"\"\n",
    "        self.p = p\n",
    "        m = n if m is None else m\n",
    "        if seed is not None:\n",
    "            np.random.seed(seed)\n",
    "        self.array = np.random.choice([0, 1], (n, m), p=[1-p, p])\n",
    "        \n",
    "        # fill the top row with wet cells\n",
    "        self.array[0] = 5\n",
    "\n",
    "    def step(self):\n",
    "        \"\"\"Executes one time step.\"\"\"\n",
    "        a = self.array\n",
    "        c = correlate2d(a, self.kernel, mode='same')\n",
    "        self.array[(a==1) & (c>=5)] = 5\n",
    "        \n",
    "    def num_wet(self, cols=None):\n",
    "        \"\"\"Total number of wet cells.\n",
    "        \n",
    "        cols: number of columns to select\n",
    "        \"\"\"\n",
    "        a = self.array[:, :cols]\n",
    "        return np.sum(a == 5)\n",
    "    \n",
    "    def num_porous(self, cols=None):\n",
    "        \"\"\"Total number of porous cells.\"\"\"\n",
    "        a = self.array[:, :cols]\n",
    "        return np.sum(a == 1)\n",
    "    \n",
    "    def bottom_row_wet(self):\n",
    "        \"\"\"Number of wet cells in the bottom row.\"\"\"\n",
    "        return np.sum(self.array[-1] == 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nothing special about the viewer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class PercolationViewer(Cell2DViewer):\n",
    "    \"\"\"Draws and animates a Percolation object.\"\"\"\n",
    "    cmap = plt.get_cmap('Blues')\n",
    "    options = dict(alpha=0.6,\n",
    "                   interpolation='nearest', \n",
    "                   vmin=0, vmax=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here an example that shows the first three time steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "seed = 22\n",
    "perc = Percolation(10, p=0.5, seed=seed)\n",
    "viewer = PercolationViewer(perc)\n",
    "\n",
    "thinkplot.preplot(cols=3)\n",
    "viewer.step()\n",
    "viewer.draw()\n",
    "\n",
    "thinkplot.subplot(2)\n",
    "viewer.step()\n",
    "viewer.draw()\n",
    "\n",
    "thinkplot.subplot(3)\n",
    "viewer.step()\n",
    "viewer.draw()\n",
    "\n",
    "plt.savefig('chap07-5.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`test_perc` runs a percolation model and returns `True` if water reaches the bottom row and `False` otherwise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_perc(perc):\n",
    "    \"\"\"Run a percolation model.\n",
    "    \n",
    "    Runs until water gets to the bottom row or nothing changes.\n",
    "    \n",
    "    returns: (boolean, number of steps)\n",
    "    \"\"\"\n",
    "    num_wet = perc.num_wet()\n",
    "\n",
    "    num_steps = 0\n",
    "    while True:\n",
    "        perc.step()\n",
    "        num_steps += 1\n",
    "\n",
    "        if perc.bottom_row_wet():\n",
    "            return True, num_steps\n",
    "        \n",
    "        new_num_wet = perc.num_wet()\n",
    "        if new_num_wet == num_wet:\n",
    "            return False, num_steps\n",
    "\n",
    "        num_wet = new_num_wet"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run a small example with `p=0.5` "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "perc = Percolation(10, p=0.5, seed=seed)\n",
    "flag, num_steps = test_perc(perc)\n",
    "print(flag, num_steps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the initial state:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "perc = Percolation(10, p=0.5, seed=seed)\n",
    "viewer = PercolationViewer(perc)\n",
    "anim = viewer.animate(frames=num_steps+1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's the animation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "anim"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a given `p` we can estimate the probability of a percolating cluster by running several random configurations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def estimate_prob_percolating(p=0.5, n=100, iters=100):\n",
    "    \"\"\"Estimates the probability of percolating.\n",
    "    \n",
    "    p: probability that a cell is permeable\n",
    "    n: int number of rows and columns\n",
    "    iters: number of arrays to test\n",
    "    \n",
    "    returns: float probability\n",
    "    \"\"\"\n",
    "    count = 0\n",
    "    for i in range(iters):\n",
    "        perc = Percolation(n, p=p)\n",
    "        flag, _ = test_perc(perc)\n",
    "        if flag:\n",
    "            count += 1\n",
    "        \n",
    "    return count / iters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At `p=0.55` the probability is low."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fraction = estimate_prob_percolating(p=0.55)\n",
    "print(fraction)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At `p=0.6`, the probability is close to 50%, which suggests that the critical value is nearby."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fraction = estimate_prob_percolating(p=0.6)\n",
    "print(fraction)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At `p=0.65` the probability is high."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fraction = estimate_prob_percolating(p=0.65)\n",
    "print(fraction)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can search for the critical value by random walk: if there's a percolating cluster, we decrease `p`; otherwise we increase it.\n",
    "\n",
    "The path should go to the critical point and wander around it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_critical(p=0.6, n=100, iters=100):\n",
    "    \"\"\"Estimate p_crit by random walk.\n",
    "    \n",
    "    returns: list of p that should wander around p_crit\n",
    "    \"\"\"\n",
    "    ps = [p]\n",
    "    for i in range(iters):\n",
    "        perc = Percolation(n=n, p=p)\n",
    "        flag, _ = test_perc(perc)\n",
    "        if flag:\n",
    "            p -= 0.005\n",
    "        else:\n",
    "            p += 0.005\n",
    "        ps.append(p)\n",
    "    return ps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see whether the critical value depends on the size of the grid.\n",
    "\n",
    "With `n=50`, the random walk wanders around 0.59."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%time ps = find_critical(n=50, iters=1000)\n",
    "thinkplot.plot(ps)\n",
    "np.mean(ps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Larger values of `n` don't seem to change the critical value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%time ps = find_critical(n=100, iters=200)\n",
    "thinkplot.plot(ps)\n",
    "np.mean(ps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%time ps = find_critical(n=200, iters=40)\n",
    "thinkplot.plot(ps)\n",
    "np.mean(ps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%time ps = find_critical(n=400, iters=10)\n",
    "thinkplot.plot(ps)\n",
    "np.mean(ps)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fractals"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Near the critical point, the cluster of wet cells forms a fractal.  We can see that visually in these examples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "seed = 22\n",
    "perc1 = Percolation(100, p=0.6, seed=seed)\n",
    "flag, num_steps = test_perc(perc1)\n",
    "print(flag, num_steps)\n",
    "viewer = PercolationViewer(perc1)\n",
    "viewer.draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "seed = 22\n",
    "perc2 = Percolation(200, p=0.6, seed=seed)\n",
    "flag, num_steps = test_perc(perc2)\n",
    "print(flag, num_steps)\n",
    "viewer2 = PercolationViewer(perc2)\n",
    "viewer2.draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "perc3 = Percolation(300, p=0.6, seed=seed)\n",
    "flag, num_steps = test_perc(perc3)\n",
    "print(flag, num_steps)\n",
    "viewer3 = PercolationViewer(perc3)\n",
    "viewer3.draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thinkplot.preplot(cols=3)\n",
    "viewer.draw()\n",
    "\n",
    "thinkplot.subplot(2)\n",
    "viewer2.draw()\n",
    "\n",
    "thinkplot.subplot(3)\n",
    "viewer3.draw()\n",
    "\n",
    "plt.savefig('chap07-6.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To measure fractal dimension, let's start with 1D CAs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from Cell1D import Cell1D, Cell1DViewer\n",
    "\n",
    "def draw_ca(rule, n=32):\n",
    "    ca = Cell1D(rule, n)\n",
    "    ca.start_single()\n",
    "    ca.loop(n-1)\n",
    "    viewer = Cell1DViewer(ca)\n",
    "    viewer.draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's one rule that seems clearly 1D, one that is clearly 2D, and one that we can't obviously classify."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thinkplot.preplot(cols=3)\n",
    "draw_ca(20)\n",
    "\n",
    "thinkplot.subplot(2)\n",
    "draw_ca(50)\n",
    "\n",
    "thinkplot.subplot(3)\n",
    "draw_ca(18)\n",
    "\n",
    "plt.savefig('chap07-7.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function creates a 1D CA and steps through time, counting the number of on cells after each time step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def count_cells(rule, n=500):\n",
    "    ca = Cell1D(rule, n)\n",
    "    ca.start_single()\n",
    "    \n",
    "    res = []\n",
    "    for i in range(1, n):\n",
    "        cells = np.sum(ca.array)\n",
    "        res.append((i, i**2, cells))\n",
    "        ca.step()\n",
    "        \n",
    "    return res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This function plots the results, comparing the rate of cell growth to `size` and `size**2`.\n",
    "\n",
    "And it uses linregress to estimate the slope of the line on a log-log scale."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import linregress\n",
    "\n",
    "def test_fractal(rule, ylabel=''):\n",
    "    res = count_cells(rule)\n",
    "    steps, steps2, cells = zip(*res)\n",
    "\n",
    "    thinkplot.plot(steps, steps2, label='d=2', linestyle='dashed')\n",
    "    thinkplot.plot(steps, cells, label='rule=%d' % rule)\n",
    "    thinkplot.plot(steps, steps, label='d=1', linestyle='dashed')\n",
    "\n",
    "    thinkplot.config(xscale='log', yscale='log',\n",
    "                     xlabel='time steps', ylabel=ylabel,\n",
    "                     xlim=[1, 600], loc='upper left')\n",
    "\n",
    "    for ys in [cells]:\n",
    "        params = linregress(np.log(steps), np.log(ys))\n",
    "        print(params[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The linear rule has dimension close to 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_fractal(20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The triangular rule has dimension close to 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_fractal(50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the Sierpinski triangle has fractal dimension approximately 1.57"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_fractal(18)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thinkplot.preplot(num=3, cols=3)\n",
    "test_fractal(20, ylabel='number of cells')\n",
    "\n",
    "thinkplot.subplot(2)\n",
    "thinkplot.preplot(num=3)\n",
    "test_fractal(50)\n",
    "\n",
    "thinkplot.subplot(3)\n",
    "thinkplot.preplot(num=3)\n",
    "test_fractal(18)\n",
    "\n",
    "plt.savefig('chap07-8.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mathematically, the fractal dimension is supposed to be:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.log(3) / np.log(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fractals in percolation models\n",
    "\n",
    "We can measure the fractal dimension of a percolation model by measuring how the number of wet cells scales as we increase the size of a bounding box."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function takes a percolation model that has run to completion.  It computes bounding boxes with sizes from 10 up to `n-1`, positioned in the center of the array.\n",
    "\n",
    "For each bounding box it counts the number of wet cells."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.stats import linregress\n",
    "\n",
    "def plot_perc_scaling(sizes, p):\n",
    "    res = []\n",
    "    \n",
    "    for size in sizes:\n",
    "        perc = Percolation(size, p=p)\n",
    "        flag, _ = test_perc(perc)\n",
    "        if flag:\n",
    "            num_filled = perc.num_wet() - size\n",
    "            res.append((size, size**2, num_filled))\n",
    "        \n",
    "    sizes, cells, filled = zip(*res)\n",
    "    \n",
    "    thinkplot.plot(sizes, cells, label='d=2')\n",
    "    thinkplot.plot(sizes, filled, label='filled', style='.')\n",
    "    thinkplot.plot(sizes, sizes, label='d=1')\n",
    "    \n",
    "    thinkplot.config(xlabel='size',\n",
    "                    ylabel='cell count',\n",
    "                    xscale='log', xlim=[9, 110], \n",
    "                    yscale='log', ylim=[9, 20000],\n",
    "                    loc='top left')\n",
    "    \n",
    "    for ys in [cells, filled, sizes]:\n",
    "        params = linregress(np.log(sizes), np.log(ys))\n",
    "        print(params[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we plot the number of cells versus the size of the box on a log-log scale, the slope is the fractal dimension.\n",
    "\n",
    "When `p` is near the critical point, the fractal dimension of the wet cells is near 1.85, but it varies from one run to the next."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sizes = np.arange(10, 101)\n",
    "plot_perc_scaling(sizes, p=0.59)\n",
    "\n",
    "plt.savefig('chap07-9.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** In Chapter 7 we showed that the Rule 18 CA produces a fractal.\n",
    "Can you find other rules that produce fractals?  For each one,\n",
    "estimate its fractal dimension.\n",
    "\n",
    "Note: the `Cell1D` object in `Cell1D.py` does not wrap around from the left edge to the right, which creates some artifacts at the boundaries.  You might want to use `Wrap1D`, which is a child class of `Cell1D` that wraps around.  It is also defined in `Cell1D.py`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** In 1990 Bak, Chen and Tang proposed a cellular automaton that is\n",
    "an abstract model of a forest fire.  Each cell is in one of three\n",
    "states: empty, occupied by a tree, or on fire.\n",
    "\n",
    "The rules of the CA are:\n",
    "\n",
    "* An empty cell becomes occupied with probability $p$.\n",
    "\n",
    "* A cell with a tree burns if any of its neighbors\n",
    "  is on fire.\n",
    "\n",
    "* A cell with a tree spontaneously burns, with\n",
    "  probability $f$, even if none of its neighbors is on fire.\n",
    "\n",
    "* A cell with a burning tree becomes an empty cell in the next\n",
    "  time step.\n",
    "\n",
    "Write a\n",
    "program that implements it.  You might want to inherit from `Cell2D`.\n",
    "Typical values for the parameters are\n",
    "$p=0.01$ and $f=0.001$, but you might want to experiment with other\n",
    "values.\n",
    "\n",
    "Starting from a random initial condition, run the CA until it reaches\n",
    "a steady state where the number of trees no longer increases or\n",
    "decreases consistently.  \n",
    "\n",
    "In steady state, is the geometry of the forest fractal?\n",
    "What is its fractal dimension?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
